A molecular dynamics study of chemical gelation in a patchy particle model 
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We report event-driven molecular dynamics simulations of the irreversible gelation of hard ellip- 
soids of revolution containing several associating groups, characterizing how the cluster size distri- 
bution evolves as a function of the extent of reaction, both below and above the gel point. We find 
that in a very large interval of values of the extent of reaction, parameter-free mean-field predictions 
are extremely accurate, providing evidence that in this model the Ginzburg zone near the gel point, 
where non-mean field effects are important, is very limited. We also find that the Flory's hypothesis 
for the post-gelation regime properly describes the connectivity of the clusters even if the long-time 
limit of the extent of reaction does not reach the fully reacted state. This study shows that irre- 
versibly aggregating asymmetric hard-core patchy particles may provide a close realization of the 
mean-field model, for which available theoretical predictions may help control the structure and the 
connectivity of the gel state. Besides chemical gels, the model is relevant to network-forming soft 
materials like systems with bioselective interactions, functionalized molecules and patchy colloids. 



o 



I 

-a 
c 

o 
o 



> 

On 

m 

(N 
O 
OO 

o 



X 



I. INTRODUCTION 



Irreversible polymerization is a mechanism of self- 
organization of molecules which proceeds via the forma- 
tion of covalent bonds between pairs of mutually-reactive 
groupsiii 2 -' 3 - If monomers with functionality (number / 
of reactive groups on a monomer) greater than two are 
present, branched molecules grow by reactions and con- 
vert the system from a fluid of monomers into a well 
connected cross-linked network, giving rise to a chem- 
ical gelation process. At the gel point, a persistent 
network spanning the sample first appears; the system 
is then prevented from flowing, yet not arrested on a 
mesoscopic length scale. The development of a network 
structure results, for example, from step polymerization, 
chain addition polymerization and cross-linking of poly- 
mer chainsi 4 ^ The same phenomenon is also observed in 
colloids and other soft materials when the thermodynam- 
ics and the molecular architecture favor the formation of 
a limited number of strong interactions (i.e., with at- 
traction strength much larger than the thermal energy) 
between different particles. Chemical gelation has been 
extensively studied in the past, starting from the pio- 
neering work of Flory and Stockmayer 1 '* 5 who developed 
the first mean-field description of gelation, providing ex- 
pressions for the cluster size distribution as a function 
of the extent of reaction and the critical behavior of the 
connectivity properties close to gelation. More appro- 
priate descriptions based on geometric percolation con- 
cepts have, in the late seventies, focused on the non-mean 
field character of the transition, which reveals itself near 
the gel point, extending to percolation the ideas devel- 
oped in the study of the properties of systems close to a 
second-order critical point. Several important numerical 
studies^ 10 ! 11 ' 12 ! 13 - 14 ! 15 ! 16 ! 17 ^ —most of them based 
on simulations on lattice — have focused on the critical 
behavior close to the percolation point, providing evi- 
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FIG. 1: Graphic description of the A and B particles (left) 
and snapshot of the simulated system (right) . The centers of 
the small spheres locate the bonding sites on the surface of 
the hard-core particle. 



dence of the percolative nature of the transition and accu- 
rate estimates of the percolation critical exponents. As in 
critical phenomena, a crossover from mean-field to per- 
colation behavior is expected close to the gel transition. 19 
But, how the microscopic properties of the system con- 
trol the location of the crossover (i.e., how wide is the re- 
gion where the mean-field description applies) and how 
accurate is the mean-field description far from the per- 
colation point is not completely understood. Another 
important open question regards the connectivity prop- 
erties of chemical gels well beyond percolation.— Even 
in the mean-field approximation, several possibilities for 
the post-gel solutions have been proposed, based on dif- 
ferent assumptions on the reactivity of sites located on 
the infinite cluster. 20 ! 21 Different propositions predict dif- 
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ferent cluster-size distributions above the gel point and 
a different evolution with time for the extent of reaction. 

Here we introduce a model inspired by step- 
wise polymerization of bifunctional diglycidyl-ether of 
bisphenol-A (B particles in the following) with penta- 
functional diethylenetriamine (A particles)^ To incor- 
porate excluded volume and shape effects, each type of 
molecule is represented as hard homogeneous ellipsoid 
of appropriate length, whose surface is decorated in a 
predefined geometry by / identical reactive sites per par- 
ticle (see Figure []}. In this respect, the model is also 
representative of colloidal particles functionalizcd with a 
limited number of patchy attractive sites,— where the 
selectivity of the interaction is often achieved building 
on biological specificit y 24 i 25 i 26 The off-lattice evolution 
of the system is studied via event-driven molecular dy- 
namics simulations, using a novel code which specifically 
extends to ellipsoidal particles the algorithm previously 
designed for patchy spheres.— Differently from previous 
studies, we do not focus on the critical properties close 
to the gel-point but study in detail the development of 
the irreversible gelation process and the properties of the 
cluster size distribution in the pre- and post-gelation 
regime. 

We find that the dynamic evolution of the system pro- 
duces an irreversible (chemical) gelation process whose 
connectivity properties can be described, in a very 
large window of the extent of reaction, with the Flory- 
Stockmayer (FS) predictions^^ This offers to us the 
possibility to address, in a well controlled model, the ki- 
netics of the aggregation and to evaluate the extent of 
reaction at which the breakdown of the Flory post-gel 
solution takes place. 



II. METHOD 

We study a 5:2 binary mixture composed of Na — 480 
ellipsoids of type A and Nb — 1200 ellipsoids of type B, 
for a total of N = 1680 particles. A particles are modeled 
as hard ellipsoids of revolution with axes a = b = 2a and 
c = lOtr and mass m; B particles have axes a = b = 4a 
and c = 20a, mass 3.4m. Simulations are performed at 
a fixed packing fraction = 0.3. Five (two) sites are 
rigidly anchored on the surface of the A (B) particles, 
as described in Fig. [1] Sites on A particles can only re- 
act with sites on B particles. Every time, during the 
dynamic evolution, the distance between two mutually- 
reactive sites becomes smaller than a predefined distance 
5 = 0.2a, a new bond is formed between the particles. 
To model irreversible gelation, once a bond is formed, 
it is made irreversible by switching on an infinite bar- 
rier at distance = 5 between the sites i and j in- 
volved, which prevents both the formation of new bonds 
in the same sites and the breaking of the existing one. 
Hence, the newly formed bond cannot break any longer, 
and the maximum distance between the two reacted sites 
is constrained to remain smaller than 5. Similarly, the 
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FIG. 2: Time dependence of the fraction of bonds p. Sym- 
bols: simulation results (averaged over 40 independent real- 
izations). For the chosen stoichiometry, p coincides with the 
reacted fraction of A reactive sites, i.e. the A conversion, or 
equivalently with the reacted fraction of B sites, i.e. the B 
conversion, p = 1 would indicate that all possible bonding 
sites have reacted. Time is measured in arbitrary units. Line: 
p(t) = kt/(l + kt), with the fit-parameter k fixing the time 
scale. This functional form is expected when any pair of re- 
active groups in the system is allowed to react, but loops do 
not occur in finite size clusters.— 

two reacted sites cannot form further bonds with avail- 
able unreacted sites. The composition of the system and 
the particle functionality are such that the reactive sites 
of type A and B are initially present in equal number, 
/aNa = IbNb, which in principle allows the formation 
of a fully bonded state in which all the sites have reacted. 
This offers a way to properly define the extent of reac- 
tion as the ratio p between the number of bonds present 
in a configuration and the maximum number of possible 
bonds JaN A - 

Between bond-formation events, the system propa- 
gates according to Newtonian dynamics at temperature 
T = 1.0. As in standard event-driven codes, the config- 
uration of the system is propagated from one collisional 
event to the next one. Note that temperature only con- 
trols the time scale of exploration of space, by modulating 
the average particle's velocity. An average over 40 inde- 
pendent starting configurations is performed to improve 
statistics. 



III. RESULTS 

In the starting configurations no bonds are present by 
construction. As a function of time, the fraction p of 
formed bonds — a measure of the state of advancement of 
the reaction — increases monotonically, until most of the 
particles are connected in one single cluster (Figure [2]). 
As a result, p saturates around 0.86, despite the fact 
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that an equal number of reactive sites of type A and B 
is initially present in the system. 

Flory and Stockmayer-i*& laid out the basic relations be- 
tween extent of reaction and resulting structure in step 
polymerizations, on the assumptions that all functional 
groups of a given type are equally reactive, all groups 
react independently of one another, and that ring for- 
mation does not occur in molecular species of finite size. 
Only when p exceeds a critical value p c infinitely large 
molecules can growi In this respect the FS theory de- 
scribes the gelation transition as the random percolation 
of permanent bonds on a loopless latticed The present 
model satisfies the conditions of equal and independent 
reactivity of all reactive sites. The absence of closed 
bonding loops in finite size clusters is not a priori im- 
plemented; as we will show in the following, however, 
such a condition — favored by the poor flexibility of the 
bonded particles and their elongated shape, the absence 
of an underlying lattice and the asymmetric location of 
the reactive sites — is valid in a surprisingly wide region 
of p values. 

The FS theory predicts the p dependence of the clus- 
ter size distribution in the very general case of a mixture 
of monomers bearing mutually reactive groups^ In the 
present case, the number n\ m of clusters containing I bi- 
functional particles and m pentafunctional ones can be 
written as 

n lm = N B N A p l+m -\l - pf m+2 w lm (1) 
(4m)! 

Wlm ~ (I - m + l)\(4m - I + l)\m\ 

and the number of clusters of size s is obtained by sum- 
ming over all contributions such that I + m — s, i.e., 
n s = Y,im,i+ m =s n i™- As shown in Figure [3k. on increas- 
ing p, the n s distribution becomes broader and broader 
and develops a power-law tail. The theory predicts a 
gelation transition when p c = (Ja — — 1) = 

0.5i£ Even close to p = 0.5, the FS prediction — which 
conforms to the prediction of random percolation on a 
Bethe (loopless) lattice where n s ~ s~ 2 5 at the perco- 
lation threshold — is consistent with the numerical data. 
On further increasing p (Figure [HJs) , the distribution of 
finite size clusters progressively shrinks, and only small 
clusters survive. Data show that Eq[TJ with no fitting pa- 
rameters, predicts rather well the numerical distribution 
at any extent of polymerization, both below and above 
the point where the system is expected to percolate, in- 
cluding details such that the local minimum at s = 2. 

To compare with the mean-field prediction of gelation 
at p c = 0.5, we examine the connectivity properties of 
the aggregates for each studied value of p, searching for 
the presence of clusters which are infinite under peri- 
odic boundary conditions. We find that configurations 
at p = 0.497 ± 0.008 have not yet developed a percolat- 
ing structure while configurations at p — 0.513 ± 0.007 
have. Hence, we locate the gel point at p c = 0.505±0.007, 
in close agreement with the theoretical mean-field expec- 
tations. Beyond this point, the material which belongs 




FIG. 3: Distribution of finite size clusters n s for different frac- 
tion of bonds p (a) below and (b) above percolation. Points 
are simulation data and lines are the corresponding theoret- 
ical curves according to FS. The dashed line represents the 
power law n s ~ s~ 2,5 . 



to the infinite (percolating) network N x constitutes the 
gel, while the soluble material formed by the finite clus- 
ters which remain interspersed within the giant network 
constitutes the sol. Figure [4ji shows that the fraction of 
gel Poo = Naa/N and even its partition between particles 
of type A (P A)OC = N Ai0a /N) and B {Pb.oo = N Bl00 /N) 
calculated according to the FS theory^ properly rep- 
resent the simulation results throughout the polymer- 
ization process. Indeed, the proportion of B particles 
to A particles in gel and in sol is a function of p (see 
inset). The relative amount of B particles in the sol 
(Nb jSO i/Na,soi) increases as a consequence of the prefer- 
ential transfer of the A particles (having more reactive 
sites) to the gel, in a way that the fraction p so i of sites B 
in the sol that have reacted (extent of reaction in the sol) 
differs from the total fraction p of sites B reacted (extent 
of reaction in the system). The constitution of the sol 
(Figure [31(b)) results to be the same as that of a smaller 
system made of Na, S oI particles of type A and Nb, S oI 
particles of type B reacted up to the extent p so i 

The evolution of the cluster size distribution can be 
quantified by the number (x n ) and weight average (x w ) 
cluster sizes of the sol, defined as x n — J2 S S? W J2s n s 
and x w = Es s2fl s/E s sn s- The numerical results and 
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fraction of bonds p 

FIG. 4: (a) Gel fraction P^ and its partition between particles 
of type A (Pa, 00) and B (Pb, 00) vs the fraction of bonds p 
(i.e. the extent of reaction in the system) . The inset shows the 
proportion of B particles to A particles in gel (Nb,oo/Na,oo 
— left axis) and in sol {Nb,soi/Na,soI — right axis) vs p. (b) 
Number and weight average cluster size (x n and x m ) prior to 
gelation and for the sol after gelation vs the fraction of bonds 
p. (c) Relation between the number of finite size clusters 
(molecules in the sol) n so i and the fraction of bonds p. The 
inset shows the number of loops ni oop vs p. In all panels, 
symbols are simulation results and solid lines FS predictions. 



the FS theoretical predictions are shown in Figure SJd. 
Both averages increase before gelation; then, they regress 
in the sol existing beyond the gel point, since large clus- 
ters are preferentially incorporated into the gel network. 
While x n increases only slightly up to the gel point, never 
exceeding 3.5, x w increases sharply in proximity of p c 
as well as sharply decreases beyond this point, consis- 
tently with the fact that x w is singular at percolation 
being dominated by large clusters. Again, simulation 
data agree very well with FS predictions. Discrepan- 
cies between theory and simulation — which reveal the 
mean-field character of the FS theory — only concern the 
range of p very near p c , suggesting that for this model 
the crossover from mean-field to percolation is very close 
to the gel point — i.e., the Ginzburg zonei^ near the 
gel point, where non-mean field effects are important, is 
very limited. A finite-size study very close to the critical 
point would be requested to accurately locate the per- 



colation point and the critical exponents, a calculation 
beyond the scope of the present work. 

From a physical point of view, the change from mean- 
field to percolation universality class is rooted in the pres- 
ence of bonding loops in the clusters of finite size, which 
pre-empts the possibility to predict the cluster size distri- 
bution. The realistic estimate of the percolation thresh- 
old and the agreement between theory and simulation 
(Fig. [3]) suggest that the present model strongly disfa- 
vors the formation of loops in finite clusters, at least 
for cluster sizes probed in this finite-size system. As a 
test, we evaluate the total number of finite (sol) clusters 
function of the extent of reaction. If 
finite clusters do not contain closed loops, n so i equals the 
number of particles in the sol minus the number of bonds, 
since each added bond decreases the number of clusters 
by one. This applies equally to the system preceding gela- 
tion, or to the sol existing beyond the gel point. Thus, at 
P < Pc (pre-gelation) the relation between n so i and p is 
linear, i.e. n so i — N — 2NbP- At p > p c (post-gelation), 
n sot can be calculated as n so i = N so i-2N B , S olPsoh where 
N so i is the number of particles in the sol fraction (Nb, so i 
of which bear reactive sites of type B) , and p ao i ^ p is the 
reacted fraction of sites B in the sol. Hence, the relation 
between n so i and p crosses to a nonlinear behavior, so 
that the number of clusters becomes one when p = 1. As 
shown in FigureUJ:, the number of finite clusters found in 
the simulation data conforms to the theoretical expecta- 
tion for all p values, both below and above the gel point. 
Hence, as a first approximation, loops are only present in 
the infinite (percolating) cluster and do not significantly 
alter the distribution of the finite size clusters, both be- 
low and above percolation. The difference between n so i 
found in simulation and the value predicted by the FS 
theory counts the number of loops in the sol, ni oop . Such 
a quantity is shown in the inset of Figure 2b. The maxi- 
mum value of nioop, achieved for p ~ p c , corresponds to 
0.2% of the total number of bonds. This demonstrates 
that intramolecular bonds within finite clusters can be 
neglected, consistent with the Flory hypothesis for the 
post-gelation regime 2 ^. Figure 2b also shows that the 
linear relation between n so i and p is valid also after the 
gel point (up to p ~ 0.6). This finding is in full agreement 
with recent experimental studie o 22 ' 31 ' 32 on the polymer- 
ization of bifunctional diglycidyl-ether of bisphenol-A 
with pentafunctional diethylcnetriamine, also suggesting 
that the number of cyclic connections in the infinite clus- 
ter is negligible well above p c . 

As a further confirmation of the absence of closed loops 
we compare the time evolution of p with the prediction 
of the mean-field kinetic modeling of polymerization, 
based on the solution of the Smoluchowski coagulation 
equation. 33 ' 34 For loopless aggregation, p(t) is predicted 
to follow 



where the fit-parameter k, which has the meaning of a 
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bond kinetic constant, fixes the time scale of the aggre- 
gation process. The time evolution of p is found to per- 
fectly agree with the theoretical predictions^ (see Fig- 
ure[2]) up to p « 0.6, i.e. beyond p c . While the prediction 
would suggest that p(t — > oo) = 1 (dash line in Figure[2]), 
the simulation shows that the formation of a percolating 
structure prevents the possibility of completing the chem- 
ical reaction, leaving a finite number of unreacted sites 
frozen in the structure. As shown above (Figure [3]), even 
in this frozen state the cluster size distribution is pro- 
vided by the Flory's post-gel hypothesis. Such a feature 
is not captured by the mean field Smoluchowski equation 
in which spatial information in the kernels are neglected. 

IV. CONCLUSIONS 

A binary mixture of patchy hard ellipsoids undergoing 
chemical gelation displays a very large interval of the ex- 
tent of reaction in which parameter-free mean-field pre- 
dictions are extremely accurate. The connectivity prop- 
erties of the model are properly described — without 
any fitting parameter — both below and above perco- 
lation by the mean-field loopless classical FS theoryiiiSI 
The mean-field cluster size distribution for the sol com- 
ponent is found to be valid for all values of the extent of 
reaction, both below and above the gel point, suggesting 
that for the present model, the Flory's hypothesis for the 
post-gelation regime properly describes the irreversible 
aggregation phenomenon, despite the explicit considera- 
tion of the excluded volume. 

The absence of loops in finite size clusters, which is 
not assumed by the model, results from the specific ge- 
ometry of the bonding pattern and by the presence of 
the excluded volume interactions, disfavoring the forma- 
tion of ordered bonding domains. Hence, the geometry 
of the particles and the location of the reactive sites on 
them may play a significant role in the stabilization of 
the mean-field universality class with respect to the per- 
colation universality class^ locating the crossover be- 
tween the two classes^ very close to the gel point. The 
present study shows that irreversibly aggregating asym- 
metric hard-core patchy particles, even if excluded vol- 
ume effects arc properly taken into account, may provide 
a close realization of the FS predictions in a wide range 
of p values. The model thus offers a starting point — for 
which theoretical predictions are available — for further 



investigations of the gelation process and for a more pre- 
cise control over the structure and connectivity of the 
gel state. In particular, a full and detailed structural in- 
formation can be known along with the dynamics of the 
system, which is potentially useful to investigate the rela- 
tion between structural heterogeneity and heterogeneous 
dynamics,— and to shed light on the microscopic aspects 
of the dynamic crossover from short 3 ^ to long relaxation 
tiniest during irreversible polymerization. 

While the structural properties are all well-described 
by the FS theory, the evolution of the extent of reaction, 
modeled via the coagulation Smoluchowski equation, is 
properly described by the theory only in the pre-gelation 
region. After gelation, kinetic constraints due to the ab- 
sence of mobility of the reactive sites anchored to the 
percolating cluster or to smaller clusters trapped inside 
the percolating matrix prevent the completion of the re- 
action and the extent of reaction freezes (to p w 0.86 
in the present case) before reaching one (as Eq. 2 would 
predict). A proper modeling of the long-time behavior 
will require the insertion of spatial information inside the 
kernels entering the Smoluchowski equation. The freez- 
ing of the extent of reaction at long times correspondingly 
freezes the cluster size distribution to that predicted by 
Flory for the reached p value. 

In the present model, the entire polymerization process 
proceeds via a sequence of FS cluster size distributions, 
determined by p(t). Recently, it has been shown that 
the FS theory properly describes also equilibrium clus- 
tering in patchy particle systems when p is a function of 
temperature and density^ It is thus tempting to spec- 
ulate that for loopless models, irreversible evolution can 
be put in correspondence with a sequence of equilibrium 
states which could be sampled in the same system for fi- 
nite values of the ratio between temperature and bonding 
depth. If this is indeed the case, chemical gelation could 
be formally described as a deep quench limit of phys- 
ical gelation. This correspondence would facilitate the 
transfer of knowledge from recent studies of equilibrium 
gel s 39 ' 40 to chemical ones. Concepts developed for irre- 
versible aggregation of colloidal particles, like diffusion- 
and reaction-limited cluster-cluster aggregation, could 
be connected to chemical gelation. Work in this direc- 
tion is ongoing. 

We acknowledge support from MIUR-PRIN. We thank 
P. Tartaglia for interesting discussions. 
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